from hep.aida import IAnalysisFactory
from hep.aida import ITreeFactory
from hep.aida import IDataPointSetFactory
from hep.aida import ITupleFactory
from hep.aida import IDataPointSet
from hep.aida import ITuple
from hep.aida import ITree
from java.lang import Math
from java.lang import String
from java.lang import Double
from java.lang import System
from java.io import File
from fr.in2p3.ipnl.aidautilities import TupleReader
from jscinttools.utilities.scintutilities import ScintUtil

thePath = ""
theFile = ""
outPath = ""

pointRT = 9
numComp = 3


treeName = "ScintillationAnalysis_ZnWO4_Txt"

af = IAnalysisFactory.create()
mainTree = af.createTreeFactory().createTree(outPath + File.separator + treeName + ".aida", "xml", ITreeFactory.CREATE)
mainTree.mkdir("Light_Yield")
mainTree.mkdir("Kinetic")
dpsf = af.createDataPointSetFactory(mainTree)
hf = af.createHistogramFactory(mainTree)
#ff = af.createFunctionFactory(mainTree)

tr = TupleReader(thePath + File.separator + theFile)
theTuple = tr.getTuple(mainTree)

temperature = []
acqWin = []
mainTree.cd("/Light_Yield")
lyDPS_Q = dpsf.create("LY", 2)
lyDPS_Corr_Q = dpsf.create("LY_Corr", 2)

decayDPS = []
lyDPS_Comp = []
lyDPS_CompRel = []

for j in range(numComp):
	mainTree.cd("/Kinetic");
	decayDPS.append(dpsf.create("Decay_Comp#" + j, 2))
	mainTree.cd("/Light_Yield");
	lyDPS_Comp.append(dpsf.create("LY_Comp#" + j, 2))
	lyDPS_CompRel.append(dpsf.create("LY_RelInt_Comp#" + j, 2))

for i in range(theTuple.rows()):
	theTuple.setRow(i);
	temperature.append(theTuple.getDouble(theTuple.findColumn("Temperature_K")))

	lyDPS.addPoint();
	lyDPS.point(i).coordinate(0).setValue(temperature[i])
	lyDPS.point(i).coordinate(1).setValue(theTuple.getDouble(theTuple.findColumn("Spectrum_mu_g")))
	lyDPS.point(i).coordinate(1).setErrorMinus(theTuple.getDouble(theTuple.findColumn("Spectrum_mu_g_Err")))
	lyDPS.point(i).coordinate(1).setErrorPlus(lyDPS.point(i).coordinate(1).errorMinus())

	amp = []
	ampErr = []
	decayT = []
	decayTErr = []

	for j in range(numComp):
		amp.append(Math.abs(theTuple.getDouble(theTuple.findColumn("Pulse_A" + (j + 1)))))
		ampErr.append(Math.abs(theTuple.getDouble(theTuple.findColumn("Pulse_A" + (j + 1) + "_Err"))))
		decayT.append(Math.abs(theTuple.getDouble(theTuple.findColumn("Pulse_t" + (j + 1) + "_us"))))
		decayTErr.append(Math.abs(theTuple.getDouble(theTuple.findColumn("Pulse_t" + (j + 1) + "_us_Err"))))
		
		if (Double.isNaN(amp[j])):
			amp[j] = 0
		if (Double.isNaN(ampErr[j])):
			ampErr[j] = 0
		if (Double.isNaN(decayT[j])):
			decayT[j] = 0
		if (Double.isNaN(decayTErr[j])):
			decayTErr[j] = 0

	acqWin.append(theTuple.getDouble(theTuple.findColumn("Pulse_Effective_window_micros")))
	corrFac = ScintUtil.calcLYCorr(amp, decayT, acqWin[i]);
	corrFacErr = ScintUtil.calcLYCorrErr(amp, ampErr, decayT, decayTErr, acqWin[i]);
	lyDPS_Corr.addPoint();
	lyDPS_Corr.point(i).coordinate(0).setValue(temperature[i]);
	lyDPS_Corr.point(i).coordinate(1).setValue(lyDPS.point(i).coordinate(1).value() * corrFac);
	lyDPS_Corr.point(i).coordinate(1).setErrorMinus(corrFac * lyDPS.point(i).coordinate(1).value() * Math.sqrt(Math.pow(lyDPS.point(i).coordinate(1).errorMinus() / lyDPS.point(i).coordinate(1).value(), 2) + Math.pow(corrFacErr / corrFac, 2)))
	lyDPS_Corr.point(i).coordinate(1).setErrorPlus(lyDPS_Corr.point(i).coordinate(1).errorMinus());

	for j in range(numComp):
		lyDPS_CompRel[j].addPoint();
		lyDPS_CompRel[j].point(i).coordinate(0).setValue(temperature[i]);

		lyDPS_Comp[j].addPoint();
		lyDPS_Comp[j].point(i).coordinate(0).setValue(temperature[i]);

		decayDPS[j].addPoint();
		decayDPS[j].point(i).coordinate(0).setValue(temperature[i]);

		if (amp[j] == 0):
			decayDPS[j].point(i).coordinate(1).setValue(Double.NaN)
			decayDPS[j].point(i).coordinate(1).setErrorMinus(Double.NaN)
			decayDPS[j].point(i).coordinate(1).setErrorPlus(Double.NaN)

			lyDPS_CompRel[j].point(i).coordinate(1).setValue(Double.NaN)
			lyDPS_CompRel[j].point(i).coordinate(1).setErrorMinus(Double.NaN)
			lyDPS_CompRel[j].point(i).coordinate(1).setErrorPlus(lyDPS_CompRel[j].point(i).coordinate(1).errorMinus())

			lyDPS_Comp[j].point(i).coordinate(1).setValue(Double.NaN)
			lyDPS_Comp[j].point(i).coordinate(1).setErrorMinus(Double.NaN)
			lyDPS_Comp[j].point(i).coordinate(1).setErrorPlus(lyDPS_Comp[j].point(i).coordinate(1).errorMinus())
		else:
			decayDPS[j].point(i).coordinate(1).setValue(decayT[j])
			decayDPS[j].point(i).coordinate(1).setErrorMinus(decayTErr[j])
			decayDPS[j].point(i).coordinate(1).setErrorPlus(decayDPS[j].point(i).coordinate(1).errorMinus())

			lyDPS_CompRel[j].point(i).coordinate(1).setValue(ScintUtil.calcLYFrac(j, amp, decayT))
			lyDPS_CompRel[j].point(i).coordinate(1).setErrorMinus(ScintUtil.calcLYFracErr(j, amp, ampErr, decayT, decayTErr))
			lyDPS_CompRel[j].point(i).coordinate(1).setErrorPlus(lyDPS_CompRel[j].point(i).coordinate(1).errorMinus())

			lyDPS_Comp[j].point(i).coordinate(1).setValue(lyDPS_CompRel[j].point(i).coordinate(1).value() * lyDPS_Corr.point(i).coordinate(1).value())
			lyDPS_Comp[j].point(i).coordinate(1).setErrorMinus(lyDPS_Comp[j].point(i).coordinate(1).value() * Math.sqrt(Math.pow(lyDPS_CompRel[j].point(i).coordinate(1).errorMinus() / lyDPS_CompRel[j].point(i).coordinate(1).value(), 2) + Math.pow(lyDPS_Corr.point(i).coordinate(1).errorMinus() / lyDPS_Corr.point(i).coordinate(1).value(), 2)))
			lyDPS_Comp[j].point(i).coordinate(1).setErrorPlus(lyDPS_Comp[j].point(i).coordinate(1).errorMinus())

for j in range(numComp):
	lyDPS_Comp[j] = ScintUtil.normalizeToRoomT(lyDPS_Comp[j], lyDPS_Corr.point(pointRT).coordinate(1).value(), lyDPS_Corr.point(pointRT).coordinate(1).errorMinus())

lyDPS = ScintUtil.normalizeToRoomT(lyDPS, lyDPS.point(pointRT).coordinate(1).value(), lyDPS.point(pointRT).coordinate(1).errorMinus())
lyDPS_Corr = ScintUtil.normalizeToRoomT(lyDPS_Corr, lyDPS_Corr.point(pointRT).coordinate(1).value(), lyDPS_Corr.point(pointRT).coordinate(1).errorMinus())

#mainTree.cd("/Kinetic");
#kineticFunction = ff.cloneFunction("Fit Function", ScintDataFitter.doKineticFit(decayDPS[2], kinParam).fittedFunction());

mainTree.commit();
